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NUMERICAL SIMULATION OF THE COMPRESSIBLE ORSZAG-TANG 
VORTEX H. SUPERSONIC FLOW 


I. INTRODUCTION 

In a recent paper, to be designated as I [1], we initiated a systematic research program 
using numerical simulations of dissipative, fully compressible magnetofiuids to identify and 
analyze the effects of compressibility on magnetohydrodynamic (MHD) turbulence. The 
first studies, reported therein, tracked the decay of initially subsonic flow fields and used 
the analogous incompressible flows as baselines for comparison. The initial average Mach 
number M ranged froi. j.2 to 0.6, and the ratio of thermal to magnetic pressure (/?) ranged 
from 30.0 to 10/3. The initial, dimensionless mass density was uniform at p = 1 while the 
Lundquist numbers varied from 50 to 200. 

The present paper deals with the effect cf embedded supersonic flows an 1 .he re- 
sulting emerging shock waves on phenomena associated with MHD turbulence, including 
reconnection, the formation of current sheets and vortex structures, and the evolution of 
spatial and temporal correlations among physical variables. The parameters used here are 
M = 1.0 and 1.5, p — 2.78 and 6.25, respectively, and /3 = 10/3, with the same range of 
Lundquist numbers as above. Often flows within regions bounded by shocks are locally 
subsonic, and the previous subsonic calculations should apply to those regions. How- 
ever, both the supersonic disturbances responsible for forming the shocks and the current 
sheets accompanying the shocks will have effects on the overall flow field which cannot 
be predicted on the basis of the subsonic studies alone. An example is the production 
of turbulence by inhomogeneous laser pulses propagating in a gas [2, 3]. The laser pulse 
produces a channel in the background gas, with some regions of the channel cross section 
heated more than ethers. The resulting shocks interact with the inhomogeneous density 
distribution of the channel to produce turbulence. This paper provides initial descriptions 
and assessments of similar supersonic effects upon a compressible magnetofluid. Our re- 
sults should be useful in describing local supersonic regions in a turbulent, compressible, 
dissipative magnetofluid. 

Among the obvious complicating factors in fully compressible MHD flows are density 
and entropy fluctuations, finite magnetoacoustic waves, discontinuities (e.g., shocks and 
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contact surfaces), and a finite information transport velocity. Theories of “weakly com- 
pressible” MHD flows [4, 5, 6, 7], while both useful and elegant, account only for density 
fluctuations (and not entropy fluctuations), filtering acoustic waves out of the solution or 
damping them out on a rapid time scale. To reveal the full range of compressible, turbulent 
phenomena, initial studies must necessarily survey the phenomenology of relevant nonisen- 
tropic flows which include all of the effects mentioned at the beginning of this paragraph. 
One can then begin to identify problems requiring theoretical treatment and exhaustive 
numerical study using the maximum achievable resolution. 

Our approach differs from that of many other studies of turbulence in that we have 
chosen a model problem, the Orszag-Tang vortex system [8], which involves decay from 
nonrandom initial conditions, rather than to follow the decay of a random initial field or 
the evolution of a driven turbuUm flow. The Orszag-Tang system is doubly periodic, and 
the initial conditions consist of single-mode, solenoidal velocity and magnetic fields, each 
containing X-points and O-points. The initial mass density is flat and the initial pressure 
fluctuations are incompressible, balancing the local forces for a magnetofluid of unit mass 
density. The identification and tracking of the various compressible MHD phenomena 
thus axe more straightforward them in a statistical approach. In addition, such a simple, 
easily reproducible, model problem can serve as a useful reference for comparison with 
the results from studies of more realistic turbulent fields [9]. This model problem, in 
fact, contains most of the significant features of MHD turbulence, including dissipation of 
magnetic and kinetic energy, reconnection, formation of high-density jets, selective decay, 
dynamic alignment, and the emergence and manifestations of small-scale structure. The 
incompressible version has received thorough study [10] and provides an excellent baseline 
for direct quantitative and structural comparisons with the evolution in a compressible 
medium. Past treatments also have emphasized the relationship of the Orszag-Tang vortex 
system to the theory of resistive tearing modes [11]. Similar highly resolved compressible 
calculations should be performed for various magnetoacoustic Mach numbers and values 
of plasma The popularity and simplicity of the model also make it an obvious choice 
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for comparison and calibration of numerical algorithms for simulating compressible MHD 
turbulence. 

Our previous paper [1] has added subsonic flow as another useful baseline for studies 
of supersonic flows. That paper provides a background on the relevant literature and a 
detailed discussion of various formulations and numerical methods used for compressible 
MHD turbulence and associated phenomena. The present paper deals with situations in 
which regions of appreciable supersonic flow lead to the emergence of shocks. These flows 
have a major effect on reconnection at the X-points in the magnetic field, the density jets 
originally observed at subsonic Mach numbers, the formation of other density extrema, 
the dissipation of energy through the formation of current sheets at the shocks, and the 
spectral composition of autocorrelations and cross- correlations of physical variables. The 
next section describes briefly our formulation and numerical methods. Section III presents 
our results on the evolution of the local structure of the flow field, the global properties 
of the system, and spectral correlations. Section IV discusses the important dynamical 
properties and observational consequences of embedded supersonic regions and emerging 
shocks in the Orszag-Tang model of an MHD system undergoing reconnection. We also 
draw conclusions regarding the effects of local supersonic regions on MHD turbulence. 


3 


II. NUMERICAL SIMULATIONS 
A. Formulation 

The following nonlinear partial differential equations, written in dimensionless form, 
govern the behavior of a two-dimensioned, compressible, dissipative magnetofluid in our 
numerical simulations: 


% = -V(pv), 

^ (pV) V ■ [pw - BB + i(p + |B| 2 )I - J-t] , 

f “ Vx[vxB-ivxB], 

^ = — V-[(E + p)v + (|B| 2 I — 2BB) • v 

_ J-V • T + J-(B • VB - VB • B) 


(la) 

(16) 

(lc) 


(Id) 


VB = 0, 


with an equation of state, 


p = ( T - 1)Z7 . 


(le) 


( 2 ) 


Here p is the mass density, v is the flow velocity, p is the mechanical pressure, I is the unit 
dyad, B is the magnetic induction field, U is the internal energy density, E is the total 
energy density, given in dimensionless units by 

£(*,() = p|v| 2 + |B| 2 + V, (3) 
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T = p/p is the dimensionless temperature, t = [(c? ; Vj + diVj) — |V • v6<j] e, ej is the 
stress tensor, and 7 = 5/3 is the ratio of specific heats. The thermal conductivity (/e), 
magnetic resistivity (7), and viscosity (/i) are constant and uniform, and we have assumed 
that the bulk viscosity is zero [12]. Important dimensionless numbers are the viscous 
Lundquist number S v = poV^Lo/p [13], the resistive Lundquist number S r = V'JjLo/ty, 
and the Prandtl number Pr = c p p/n. 

In these definitions, po is a characteristic density, VJJ is a characteristic Alfven speed, 
L 0 is a characteristic length equal to the reciprocal of the minimum wavenumber fc m i n = 
27 r/X (L = system length), c p is the specific heat at constant pressure, and Vo = VjJ is a 
characteristic flow speed. The thermodynamic normalization sets 

E 0 =po=Bl/Sn = p 0 vf/2 . (4) 

Time (<) is measured in units of the Alfven transit time Lq/V%. For the runs reported 
here, we set S v = S r = 50, 100, or 200 and Pr = 1 All of our numerical simulations 
assume unit magnetic Prandtl number, t.e., p/(poij) = 1, so we subsequently consider only 
the parameter S = S v — S r . 

B. Initial Conditions 

The initial magnetic and velocity fields are identical to the vortex system of Orszag 
and Tang (1979): 


and 


p(x,y,t = 0) = constant = p° , (6a) 


v(i,y,t = 0) = — siny e* +sinx e„, (65) 


B(z, y,t = 0) = — siny e* + sin2x e v , (6c) 
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where e* and e y are unit vectors in the x and y directions respectively. For reasons given 
below, our investigations of emerging shocks use an initial mass density which can differ 
from unity. The initial cross helicity and normalized correlation between the velocity and 
magnetic fields then equal 0.5. As in our previous calculations, the initial mechanical 
pressure, p(t = 0), comprises a mean part, (p(t = 0)) s p‘ , and a fluctuating part, 
6p(x, y, t = 0) = t>p°(x, y ). Here the bracket notation represents a volume average, i.e., an 
average of a function (/) over the entire grid (Q): 

J f(x,y,i) dx dy 

w - • (7) 

a 

The equation for the dimensionless speed of sound, 

CS = g , (8) 

constrains the toted pressure to equal or exceed zero. The lower bound on p° is then a 
positive definite functional of Sp°(x,y), given by p£, in > max{— 6p 0 }, assuming that the 
choice of 6p° is not identically zero in the computational domain. 

Equation (8) indicates that formula for the local Mach number M is 


= ^ . 

7 P 

Notice that the pressure must be greater than zero everywhere in order for M to be 
real and free of singularities throughout the grid. To obtain a “characteristic” or an 
“average” Mach number M at time t = 0 for the flows in our calculations, we substitute 
p = p°, p = p°, and |v°| = y/{v 2 (t = 0)) into Eq.(9a): 

^ = . 

7 P° 

The ratio of thermal to magnetic pressure depends on local Mach number according to 
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( 10 ) 


q~ P _ 2 P\ y \ 2 
P ~ |B| 2 7^ 2 |B|2 • 

In all calculations, the initial “characteristic” value, /? c , is 10/3, to permit comparison with 
the M = 0.6 case considered previously [1]. The quantity 0 C relates to 0 in Eq.(10) as M 
does to M. 

Equation (9a) shows that one could obtain flows with supersonic regions merely by 
setting the initial background pressure p° close enough in magnitude to the minimum 
(negative) value of the fluctuating pressure 6p°. As in our subsonic calculations, 6p° 
would be the solution of the incompressible MHD equations with p° = 1. Unfortunately as 
p° — ► max{— 6p 0 }, the supersonic regions move inward toward the center of the grid, which 
is a stagnation point of the flow and an X-point of the magnetic field. Figure 1(a) shows 
the initial local Mach number JA when the minimum value of the thermal pressure is 10~ 4 . 
The average background pressure is p° « 2.6, making M « 0.679. The initial maximum 
value of M is approximately 2.2, since the minimum pressure occurs near the stagnation 
point. Calculation of this case up to t = 8 shows only small differences from the M — 0.6 
subsonic case studied previously, even though the latter had an initial maximum value of 
0.97 for M. In fact, both calculations exhibit supersonic flows early, but no shocks occur 
[1]. Erlebacher et al. have reported a similar result for 3-D, compressible, Navier-Stokes 
turbulence [14]. 

For the study of effects related to emerging shocks within this model, another prescrip- 
tion is necessary. Our approach is to increase p° above 1.0 while retaining the definitions of 
v° and B° in Eqs.(6b) and (6c). In addition, we allow the fluctuating pressure, £p°, to sat- 
isfy the incompressible equation of motion with p° set equal to one. With the proper choice 
of p°, we can define a set of initial conditions which are the same as those of the previous 
subsonic study except for the initial mass density. By increasing p°, we can raise the char- 
acteristic Mach number (Eq.(9b)) of the flow above 1.0. Since these choices of p° and 6p° 
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do not satisfy the incompressible equation of motion, an imbalance exists between the mo- 
mentum gradients of the fluid and the opposing magnetic and thermal pressure fields. This 
imbalance, in fact, causes the supersonic disturbances which we desire to study. In support 
of this approach, Appendix A shows that attempting to maintain a balance among local 
forces confines us to subsonic characteristic Mach numbers M < 0.679, if the amplitudes 
of the velocity and magnetic field fluctuations are equal, as in Eq.(6). 

Our chc of initial average pressure p° is 10/3, and the corresponding characteristic 
0c is also 10/3. These values are the same as in the M = 0.6 subsonic calculation [1], in 
which p° = 1. The present calculations, therefore, are linked directly to our previous work. 
To ensure that significant supersonic regions are present, we set the characteristic Mach 
number M > 1, giving an initial ma s density p° > 2.78 by Eq.(9b). Here we calculate 
the system evolution for M = 1.0, p° = 2.78, and M = 1.5, p° = 6.25. Figure 1(b) shows 
the initial form of the local Mach number, At, for the M — 1.0 calculation. The local 
Mach number varies widely over the system from a value of 0.0 at the stagnation point 
in the center of the computational grid to a value of 1.6 within sizable supersonic regions. 
Section III shows that these regions produce shocks which alter the evolution significantly 
from the subsonic and incompressible baselines. 

C. Numerical Algorithm 

Previous papers have described our numerical algorithm in detail [1, 15). Our code, 
CRUNCH2D, implements a Fourier collocation method and employs an isotropic trunca- 
tion in Fourier space at each time-level [14]. The modified Euler method, a second-order 
Runge-Kutta scheme [16], discretizes time. The time step, At, is limited by a compressible 
MHD Courant-Friedrichs-Lewy (CFL) number, 

Mcfl = sup{ (|v| ^ } < 0.3 , (11) 

where Cs = sound speed and h = 2x/N = Ax = Ay. When shocks are present, reduction 
in the time step according to Eq.(ll) (as compared to subsonic flows) often is insufficient 
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to guarantee a smooth solution in the vicinity of the shock. In such situations, the thermal 
pressure and the mass density can have local negative values, and increased viscous and re- 
sirtive dissipation of^-n are necessary. When the diffusive terms are considered separately, 
numerical instability can occur due to enhanced transport in the vicinity of the shock 
discontinuity. We have, therefore, added stability criteria based on viscous dissipation [17| 
and resistive dissipation, respectively: 


M v = 


Mr 


N*At 

Pmia 


< 1 


45r “ 


1 


( 12 ) 


Our simulations consume approximately 17 ps per tune step per grid point on a single 
processor of the Numerical Aerodynamic Simulation Program (NAS) Cray Y-MP computer 
at the NASA Ames Research Center. For M = 1 and Lundquist numbers S = 50, we use 
256 2 collocation points while for M -- 1, S > 100, and A# = 1.5, S — 50, we use 513? 
collocation points. The calculations with 512 2 collocation points r eq uire approximately 7 
million words of core memory. Wi th these dissipation levels and resolutions, CRUNCE2D 
can simulate processes with moderate strength diodes (M < 3.0), as verified by a member of 
standard tests, e.g., the rupturing diaphragm problem [15] and the numerical anwhlioos in 
I. Our model shocks are approximately 5 computational cells thick when S = 200 and N = 
256 2 [15]. 
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ITL RESULTS AND DISCUSSION 


As in oar previous paper fl], our investigation of compressibility in the Orazag-Tfcng 
vortex s»du refies primarily on such diagnostics as time-dependent area integrals, contour 
diagrams, gesacraimed streamline plots. »nd wavenumber spectra of various autocorrelations 
and cross oo'-dataons. The render should refer to that paper for the correspo n ding {Jots 
when V r an ges foam 0-2 to 0.6. Supersonic compressioa and the presence of shocks pro- 
foundly jafiafa ce energy dissipation, the time at which the onset of strong reconnection 
occurs at the magnetic X-poaut, the longevity of reconnection and the associated current 
sheet, and the ecobUot and extent of the correlation between the magnetic and velocity 
adds. The w a ve num ber spectra and the evolution of global aver a ge s thus differ signifi- 
cantly focan those of subsonic floors. To aid in the interpretation of those differe n ce s , the 
dxacasson begins with structural features. 

A. Structural Effects of Supersonic Flows and Emerging Shocks 

For tbe in comp ressible and subsonic sys tems, the structures of interest are those as- 
sociated with magnetic field reconnection: magnetic islands and their coalescence, jets 
emerging from the reconnection zone, vortex quadrupoles, and electric current sheets (“re- 
conseciioo current sheets'" or RCS). 1c the supersonic cases discussed below, shocks and 
their associated current sheets (‘‘"shock current sheets" or SCS) and distributions of vortic- 
ity and dilatation (V - v) constitute additional structural features which become dominant 
as Mach number increases. Since the system is decaying with time, the most interesting 
structural effects occur for times t < 4.0. Most of the results presented in this subsection, 
therefore, correspond to that temporal range. 

An additional process emer g e s due to the self-interaction of the magnetofiuid via shock 
wave propagation through the system. Because of the symmetry and periodicity of the 
Orszag-Taag model, the shocks (and later, finite acoustic waves) can interact multiple 
times; with each other and with the spatially varying mass density field. The latter inter- 
actions generate progressively smaller-scale vortex structures and density fluctuations, and 
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can cause the magnetofluid to appear turbulent at later times [1, 18, 19]. Picone tt aL [20] 
have provided a simple numerical demonstration of this phenomenon. The newly generated 
vorticity affects the local mass density on time scales corresponding to the induced rota- 
tional motion rather than the shock transit. Thus the readily detectable manifestations 
appear gradually and with temporally increasing prominence. 

To illustrate the shock-emergence process, we follow the time development of the lo- 
cal Mach number M throughout the grid. Figure 2 shows the evolution of M for the 
calculation in which the average initial Mach number M = 1.0 and 5 = 100. The four 
regions of supersonic flow visible in Fig. 1(b) have moved clockwise and outward, com- 
pressing the fluid ahead of them and forming shocks around t = 1.5. Because of the 
periodic boundaries, identical flows squeeze the fluid from opposite directions, producing 
broad regions of compression. The local Mach number peaks at around 2.0 and 3.0 for 
M = 1.0, S = 100, and M - - 1.5, 5 = 50, respectively. The shock boundaries expand, 
eventually intersecting at the central X-point (t = 2.5), compressing the weak reconnec- 
tion current structure residing there (Fig. 3), and producing complex transmitted shock 
structures (t > 2.5). The transmitted shocks propagate toward the boundaries, interacting 
with the “parent” shocks and weakening through dissipation and expansion. 

The emerging shocks strongly influence the spatial variation of the magnetic field. 
Figure 3 shows the evolution of the electric current density over the interval t = 2 to t = 
4 for M = 1.0, S = 100. Current sheets (SCS) are coincident with all shocks; the 
sharpest SCS reside at the two strongest shocks. Unlike the subsonic and incompressible 
sys tems , therefore, significant dissipation of magnetic energy occurs at regions other than 
the vicinities of magnetic X-points. The magnetic energy available to drive reconnection at 
the X-points is thus lower than for M < 0.6. At t = 2, a wide, weak reconnection current 
sheet (RCS) is present in the center of the grid, where strong reconnection and sharp RCS 
were observed much earlier in the subsonic and incompressible simulations. Compression 
of this central region occurs when the shocks intersect at < = 2.5. A sharp RCS appears at 
that time, indicating that significant reconnection and resistive heating are taking place. 


11 


The current density peaks at t « 3, later than in the subsonic cases, and the central RCS 
is prominent through t — 6. (Note that weak current sheets also appear in the vicinity 
of the shocks transmitted through the center of the grid.) In the M = 1.5 calculation, 
reconnection is delayed more significantly, reaching a peak at approximately t = 4. 

The dilatation evolves similarly to the shock current sheets while significant vortex 
structures reside with both the reconnection current sheet and the shocks. Figure 4 shows 
the time development of the dilatation. High negative values near the shocks indicate that 
strong compression is occurring there. The weakening compression with time, evident in 
the decreasing minimum values, is due to the relatively low Lundquist numbers of the 
simulations and to the expansion of the shocks that were transmitted through the central 
reconnection region. Vortex structures occur in regions of greatest shock curvature, and 
an additional quadrupolar vortex resides at the central RCS, as in the incompressible 
and subsonic calculations [1]. The quadrupolar vortex structures experience increasing 
distortion with increasing Mach number. Figure 5 shows that, when the shocks reach the 
center of the grid at t = 2.5, shock vortex structures diagonally bisect the reconnection 
vortex quadrupoles. 

In the incompressible and subsonic systems, the alignment (or correlation) of the ve- 
locity and magnetic fields grows over time from the initial value of 50% (cf. I). The onset 
of this growth occurs slightly later in the subsonic simulations, the delay increasing with 
M. The next subsection shows that the supersonic cases undergo similar growth of cor- 
relation after much larger initial delays. The period of growth ends with a plateau or a 
period of much slower growth. The correlation value at the plateau is lower as Mach num- 
ber increases. Retardation of reconnection, as demonstrated by the above contour plots, 
appears to be a significant factor in reducing the rate of growth and level of correlation 
between the velocity and magnetic fields relative to the cases for which M < 0.6. Figures 
6(a) and (b) show streamline and magnetic field line plots at t = 3 for M = 1.0, S = 50. 
The magnetic field has been compressed along the horizontal direction and stretched in 
the vertical direction relative to the subsonic systems. The shocks themselves cause a kink 
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in the magnetic field, but their presence is much less apparent than in the electric current 
density, which depends on field gradients. The field lines are bent toward the shock fronts; 
hence, we are observing “fast” shocks. The velocity field at t = 3 shows the early stages of 
reconnection similar to the situation at t = 1 for M < 0.6. In particular, an X-point forms 
in the velocity field near the center of the grid, where the X-point in the magnetic field 
resided at time t = 0. This X-point in the velocity field remains in the center of the grid 
until t = 6. By that time, dissipation has significantly reduced both the intensity of recon- 
nection and the effects related to shocks. In Figs. 6(c) and (d), the velocity and magnetic 
fields at t — 7 exhibit mainly a large-scale similarity. We address this more thoroughly in 
the next subsection through a discussion of the correlation between the momentum and 
magnetic fields. 

The explanation for the delay and suppression of reconnection lies partly in the out- 
ward movement of the supersonic regions due to the initial imbalance between the cen- 
trifugal force of the magnetofluid and the opposing thermal and magnetic pressures. This 
rarefies the region near the central magnetic X-point and opposes the flow of fluid into 
the reconnection region. The time at which shocks reach the center of the grid marks 
the arrival of appreciable fluid from the outer regions. Only then does one observe rapid 
reconnection and energy dissipation, as seen at earlier times for M < 0.6. The presence of 
current sheets at the shocks also weakens reconnection by rapidly dissipating magnetic en- 
ergy that would otherwise be dissipated at the magnetic X-point. The transmitted shocks 
propagating away from the reconnection region, although weak, impede the flow of fluid 
into that region, somewhat reducing the intensity of reconnection. 

In summary, our simulations show that the emergence of shocks outside of a magnetic 
X-point can impede reconnection at the X-point through Ohmic dissipation of magnetic 
energy and through rarefaction of the reconnection region itself. Shocks reaching the 
X-point compress the reconnection region and signal the arrival of higher density mag- 
netofluid, accelerating the reconnection process. We expect the emergence of compressed 
regions bounded by shocks to affect reconnection similarly in a realistic turbulent MHD 
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flow with local supersonic regions or pressure-momentum imbalances. In a more general 
situation, reconnection might be enhanced in one region by compression while being de- 
layed by rarefaction or shock-related Ohmic dissipation in a neighboring area. 

Figure 7 shows the evolution of the mass density. Note that the system exhibits 180° 
rotational symmetry, as does the initial magnetic field. Fluid and magnetic field squeezed 
by the four initial supersonic regions (Figs. 1 and 2(a)) have been pushed rapidly along a 
direction tangential to the bounding shocks (those which emerged at t = 1.5) and inward 
toward the middle of the system. At t = 2, the shocks at the upper left-hand (ULH) 
and lower right-hand (LRH) regions extend to the central region, and the remaining shock 
structures appear to include Mach stems. The magnetofluid along that diagonal (ULH 
shock - central region - LRH shock) is rarefied while the compression regions in the other 
quadrants consist of two narrow plateaus separated by a narrow region of somewhat lower 
density. At t = 2.5, the ULH and LRH shocks have intersected at the center, compressing 
the reconnection current sheet and the magnetofiuid. Th«, circular patterns at the left 
and right ends of the RCS correspond to the ends of the fluid jet. In Fig. 7(c) (t = 3), 
the comers of the transmitted shocks and intersections of the primary or incident shocks 
have introduced additional local density peaks in the region of the jet. By t = 4, the 
density extrema have coalesced and spread and the transmitted shocks have expanded and 
weakened, interacting with the primary shocks that preceded reconnection. Figure 7(e) 
shows the mass density at the end of the calculation (t = 8). The maximum and minimum 
values indicate that the mass density is evolving toward a uniform state, consistent with 
the dissipative nature of the medium. However, because of multiple shock interactions with 
the density field during the system evolution [20] and because of the injection of small scale 
fluctuations by magnetic reconnection, the fluid has a turbulent appearance with sizable 
fluctuations at a range of scales, as mentioned at the beginning of this section. 
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B. Evolution of Global Properties 
B.l. Mach Number Dependence 

The M = 0.6, 5 = 50, calculation described by Dahlburg and Picone [1] provides a 
useful baseline for evaluating the evolution of global properties related to dynamic align- 
ment, to production of vorticity and electric currents, and to energy transfer among MHD 
variables due to reconnection, dynamo action, and shocks. Figure 8 shows the evolution 
of several global parameters for the simulations with S = 50 and M = 0.6, 1.0, and 1.5. 

Figure 8(a) shows the compressible coefficient of correlation between the momentum 
and magnetic fields: 


a c 


(pv ' B) 

V((p|v]) 2 ") (|BP> ‘ 


(13) 


For the compressible, supersonic, Orszag-Tang vortex system studied here, a c = 0.5 at 
t = 0. As indicated previously [1], the growth of a c relates to the process of “dynamic 
alignment” in turbulent, compressible magnetofluids. In incompressible models, dynamic 
alignment occurs as the kinetic and magnetic energy decay selectively with respect to the 
cross helicity. The measure often used in incompressible calculations is, therefore, the ratio 
of the cross helicity ( H e ) to incompressible total energy (Ei), given by 


,_Hc 2(v • B) 

€ Ei (|v|*) + {|B|*) • 

In all of our investigations, the two measures of correlation have qualitatively similar 
evolutions and provide useful measures of differences among incompressible, subsonic, and 
supersonic MHD flows. 

Figure 8(a) shows that the baseline subsonic system (M = 0.6) undergoes more rapid 
growth in a c than do the supersonic cases ( M > 1)). In fact, the M = 0.6 calculation has 
two stages of rapid growth (for t < 8), beginning at approximately t = 1.5 and 4.5, respec- 
tively, and ending with plateaus at successively higher correlation values. The correlation 
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q c in each of the supersonic calculations exhibits only a single stage of rapid growth to a 
plateau for t < 8, beginning at t zz 3.5 when M = 1.0 and at (« 4.5 when M = 1.5. The 
supersonic systems also experience an early stage in which the correlation level wanders 
around the initial value of 0.5 before the onset of rapid growth. The end of this initial 
“induction period” occurs at roughly the same time as the peak in the current density in 
the reconnection region (Section III. A.). Incidentally, the M = 0.6 curve shows an initial 
period of slow growth which also ends when the current density peaks in the reconnection 
region, in this instance at < « 2.0. Figure 8(a) shows that the plateau value of a c at late 
times varies approximately linearly with the reciprocal of the characteristic Mach number 
M. 

In the compressible Orszag-Tang vortex model, growth of correlation thus appears 
to depend strongly upon the dynamics of, and energy dissipation by, reconnection at the 
central magnetic X-point. As indicated in the previous subsection, the supersonic flows 
in the Orszag-Tang vortex initially rarefy the plasma in the central region, weakening 
reconnection there. Shocks bound the regions of high density, so that the arrival of shocks 
in the reconnection region is a precursor to considerable inflow of magnetofluid. The latter 
inflow is, in turn, necessary for significant reconnection to occur, and o c begins sustained 
growth shortly thereafter. The increased delay in sustained growth of correlation with 
increasing Mach number thus corresponds to the later arrival of the shocks at the central 
X-point. In addition, as the characteristic Mach number increases, the shocks become 
stronger, as do the current sheets associated with them. Then the amount of magnetic 
and kinetic energy transformed into internal energy at the shocks also must increase with 
M. This leaves less energy to drive reconnection and, consequently, dynamic alignment, 
causing the value of a c at a given time to be lower as the Mach number increases. Following 
the rapid growth of correlation, all cases reach a stage of very slow growth, after the 
reconnection event has run its course. 

Figures 8(h) and (c) show the temporal evolution of the kinetic energy (normalized to 
onej and magnetic energy for various Mach numbers. The dependence of both on M > 1 
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is much stronger than for the subsonic flows ( M < 0.6), which showed little variation with 
Mach number. In incompressible numerical simulations [10], the oscillations in kinetic and 
magnetic energy were approximately ISO 0 out of phase and were interpreted as Alfvenic. 
Our subsonic calculations [1] were quite similar to the incompressible calculations in this 
regard, and the same interpretation applied, since the initial kinetic and magnetic energies 
were equal and smaller than the internal energy' by a factor of at least five and because 
the effects of acoustic waves were small. In the present calculations, however, the initial 
kinetic energy increases as M 2 while the initial magnetic energy remains at unity. In 
addition, shocks and magnetoacoustic waves are present, so that the temporal evolutions 
of the magnetic and kinetic energies are not linked in the same manner as for the subsonic 
and incompressible systems. The increased magnitude of the kinetic energy as M increases 
explains the greater relative sensitivity of the magnetic energy to supersonic Mach number 
(Fig. 8(c) vs. 8(b)). Small percentage changes in kinetic energy could result in larger 
relative changes in the magnetic energy, whereas small fluctuations in the magnetic energy 
would affect the kinetic energy' even less. 

Additional phenomena should affect the evolution of the magnetic energy' shown in 
Fig. 8(c). First, appreciable reconnection is delayed and apparently weakened as the Mach 
number increases. Second, given that we are observing fast shocks, compression should 
enhance the magnetic energy, leading to larger peak values with increasing Mach number. 
Because of the larger jump in tangential magnetic field across the shock with higher Mach 
number, the dissipation of magnetic energy by the shocks and their associated shock current 
sheets also should increase with Mach number. In Fig. 8(c), the peak magnetic energy 
and the subsequent dip are both proportionately larger for M — 1.5 than for M = 1.0. 
In addition, the final stage of decreasing magnetic energy occurs later because the stage 
of strong reconnection is delayed. Notice that a shoulder appears at t ~ 3.5 for M = 1.0 
and at t % 4.5 for M = 1.5. These are the respective times at which the rapid growth of 
correlation began. 
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In I, our interpretation of the temporal evolution of the M = 0.6 calculation was 
different from that in the previous paragraph because no shocks were present. In Fig. 
8(c), the magnetic energy is enhanced at t « 2 when M = 0.6. This is later than for 
the supersonic cases and is consistent with an Alfvenic exchange of energy from the mo- 
mentum field to the magnetic field, similar to that occurring in incompressible media [1, 
21]. The subsequent decay of magnetic energy is more immediate than in the supersonic 
calculations, and no “shoulder” appears in the curve. We attribute this to the fact that 
strong reconnection is not delayed when M is subsonic [1]. An Alfvenic transfer of energy 
between the momentum and magnetic fields is likely to occur in the supersonic systems as 
well, and we cannot state with certainty that compressions! effects are more important. 
However, the stronger Mach number dependence of the supersonic cases leads us to believe 
that the interpretations in the previous paragraph are reasonable. 

Figures 8(d) shows the evolution of the kinetic enstrophy, (a; 2 ), where u> is the vor- 
ticity. The decrease and delay in kinetic enstrophy production relative to the subsonic 
calculation (M = 0.6) must relate to shock formation and the delay in appreciable recon- 
nection, as discussed previously. Here the initial rise in enstrophy occurs when shocks first 
emerge, indicating that shock curvature is contributing significantly to vorticity produc- 
tion, according to Crocco’s Theorem [22]. The structure in the supersonic systems depends 
in a complicated manner upon shock formation and interactions, as well as upon vorticity 
production at the central reconnection current sheet and the secondary RCS residing on 
the boundaries. 

Figure 8(e) shows the magnetic enstrophy, (j 2 ), where j is the electric current density. 
The effects of shock emergence on magnetic energy dissipation appear prominently in 
the form of a massive enhancement at t « 1.5 for M = 1.5 and as a smaller peak for 
M = 1.0. The supersonic cases have similarly sized peaks at t « 3 ( M = 1.0) and 
t « 4 (M = 1.5), corresponding to reconnection at the central magnetic X-point. In 
contrast, the broad enhancement in the subsonic system indicates that steady reconnection 
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occurs at approximately the same level over the interval 1.0 < t < 3.0. The final stage of 
decay occurs later as the Mach number increases. 

B.2. Lundquist Number Dependence 

Figure 9 shows global averages corresponding to M = 1.0 and 5 = 50, 100, and 200 
for 0 < t < 8.0. In Fig. 9(a), the compressible correlation coefficient, o c , shows the same 
qualitative behavior as in the analogous subsonic and incompressible simulations [1], with 
the degree of alignment at a given time decreasing as S increases. The kinetic energy (Fig. 
9(b)) also displays similar behavior to its subsonic counterparts, with almost no Lundquist 
number dependence during the early stage of development. However, the timing of the 
fluctuations in average kinetic energy is different from that of the subsonic calculations 
(Fig. 5(c), I). The fact that the initial kinetic energy in Fig. 9 is nearly three times the 
initial magnetic energy means that the latter affects the kinetic energy much less than in 
the subsonic flows. Apparently the compression of the fluid downstream of the initial high 
Mach number regions and the formation of shocks reduce the kinetic energy and transfer 
energy to the magnetic field. This energy transfer occurs through Alfven waves in an 
incompressible magnetofluid. 

The evolution of magnetic energy as a function of Lundquist number (Fig. 9(c)) 
is less similar to that in the subsonic cases. The initial extremum in magnetic energy 
occurs earlier, at t « 1.5, when shocks form, as opposed to t « 2 for M < 0.6. For 
S = 100 and 200, a second extremum appears at f «4 which is not seen for 5 = 50 or for 
M < 0.6. This is most likely due to a combination of two factors: 

(1) Following the arrival of the “incident” shocks at the central magnetic X-point and the 
onset of strong reconnection there, transmitted and diffracted fast shock waves interact 
with the reconnection current sheet and with remnants of the incident waves. The array 
of shocks no longer forms a network bounding quasistationary regions, as seen at earlier 
times. Instead the shocks expand and propagate through a large portion of the system, 
compressing the magnetic field and increasing the magnetic energy. Interactions among the 
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various shocks and magnetoacoustic waves enhance the complexity of the wave structure 
over time. The coverage of the grid by these waves thus increases with time for t > 1.5, and 
their influence later diminishes as dissipative effects become more prominent for t > 4.0. 
This sequence is implicit in the broad peak observed in Fig. 9(b). 

(2) The reduced rate of dissipation of shock waves and other magnetoacoustic waves at 
higher Lundquist numbers increases the strength of the compression. The broad peak in 
magnetic energy at t fa 4 thus becomes more prominent as Lundquist number increases. 

Figures 9(d) and (e) show that the kinetic and magnetic enstrophies increase approxi- 
mately linearly with Lundquist number and have similar structures. Distinct peaks emerge 
at t « 3 and t « 4 as dissipation decreases. The peak at t = 3 corresponds to the time at 
which the current in the reconnection region reaches a maximum. As in Fig. 9(c), the peak 
in current density at t = 4 reflects the wide distribution and interactions of numerous mag- 
netoacoustic waves and shock waves, following the arrival of the “incident” shock waves 
at the central magnetic X-point. The subsequent onset of strong reconnection and the 
transmission and diffraction of the incident shocks through the reconnection region (e.g., 
Fig. 2(d)) results in a complex global magnetic-field structure with significant gradients 
throughout the system at t « 4. Figure 3(d) verifies this, showing the presence of sizable 
currents over a large portion of the grid at that time. 

C. Spectral Decompositions of Correlations 

Here we examine the calculation with M = 1.0, S = 100, although the trends iden- 
tified are consistent across all of the simulations under discussion. Prior to the formation 
of shocks at t « 1.5, the wavenumber spectra are typical of decay problems and resemble 
those of the subsonic cases [1]. Figure 10(a), for example, shows the mass density auto- 
correlation spectrum at t = 1. However, after shocks form, the wavenumber spectra of the 
mass density, velocity, nonsolenoidal and solenoidal components of the velocity, and the 
magnetic field appear quite changed, as shown in Fig. 10(b) and (d) - (g) for time t = 2.0. 
Except at wavenumber extremes, each is approximately linear, with a slope similar to that 
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of the reference line (k~*). The magnetic field deviates from the others, showing a slightly 
steeper slope. The magnetic field spectrum also differed from that of the other functions 
in the subsonic simulations. 

This change of spectrum coincides with the emergence of compressed regions bounded 
by shocks. With the passage of time, the high wavenumber portion of the spectrum falls 
from the reference line as a result of dissipation. The density spectrum in Fig. 10(c) exem- 
plifies that effect. In an effort to model the spectral features of these compressed regions, J. 
Dahlburg [23] has computed numerically the equivalent spectra of elevated plateaus with 
either square or circular cross sections. The circular enhancement produced an approx- 
imate k~ 3 au ■: correlation spectrum. The square plateau had two distinct components: 
A: -2 , characteristic of a one- dimensional planar shock, and approximately k~ 3,s . Since 
the actual structures have cross sections similar in a very coarse view to a quadrilateral, 
the square “model” is closer to the simulations. These considerations show that simple 
idealizations of the two-dimensional compressed regions are not entirely adequate to ex- 
plain our observations quantitatively. However, the qualitative results do seem to explain 
the change in autocorrelation spectra which we observe when two-dimensional compressed 
regions with bounding shocks appear. 

Figure 11(a) and (b) show the spectral amplitude of the cross-helicity normalized by 
the spectral amplitude of the incompressible energy at t = 2.0 and 6.0, respectively. In the 
subsonic cases, the low wavenumber region had positive values while the high wavenumber 
region showed an anticcrrelation, although this bifurcation was less pronounced at M = 0.6 
than at M = 0.2. Here we see that the velocity and magnetic field are approximately un- 
correlated >~' , er most of the spectrum at times prior to the interaction of the shocks near 
the X-point. For later times (t > 4) with M = 1.0, the lower and middle wavenumbers 
show a positive correlation while the higher wavenumber comp nents are approximately 
uncorrelated. For M = 1.5, 5 = 50, the middle and higher wavenumber portions of the 
spectrum are uncorrelated, while the lower wavenumbers have positive values through- 
out the calculation. The presence of two-dimensional compressed structures bounded by 
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shocks, the coincident dissipation of energy by the shock current sheets, and the delay and 
reduction of energy dissipation by reconnection are the likely causes of the reduction of 
correlation at the smaller scales. In addition, these factors probably are responsible for 
the lower global (average) correlation between the magnetic and velocity fields relative to 
subsonic flows. 

Mechanisms which suppress the growth of correlation between the velocity and mag- 
netic fields at temporal scales of a few hours or less are of interest in interpreting mea- 
surements of the solar wind [24], The above compressive mechanism can occur at the 
scalf of any nonuniformity within, or encountered by, the solar wind. This is because such 
n muniformities can generate shocks, given that the solar wind is supersonic and is not a 
steady flow. Presently the excitation of Kelvin-Helmholtz instability by velocity shear is 
the leading candidate for the cause of a reduction in correlation between the magnetic and 
velocity fields with increasing distance from the Sin. t 24]. We point out that the interpre- 
tation of the data is still uncertain. For example, early data analysis indicated that MHD 
fluctuations in the solar wind are mainly Alfvenic and propagate away from the Sun [25]. 
Subsequent theoretical work inferred from this observation that the population of either 
the aligned or the antialigned state would become dominant with increasing distance from 
the Sun (e.g., [26]). Recent data analysis has indicated a much smaller population of out- 
wardly propagating Alfvenic fluctuations near the earth. This has led to tne conclusion 
that neither state of alignment dominates at distances around 1 AU from the Sun. 

The correlation of the fluctuations in thermal and magnetic pressures is mostly neg- 
ative soon after the shocks have formed, as exemplified by Fig. 11(c). At later times 
(t > 3.0), this effect becomes less pronounced or disappears entirely, depending on M and 
S. According to arguments by Barnes [27], a negative correlation is consistent with the 
presence of quasistationary structures in the flow field. Clearly such structures exist here 
until the bounding shocks reach the central magnetic X-point, after which time the shocks 
weaken and the original compressed regions change significantly. This could explain the 
above observation that the negative correlation disappears at later times. In their studies 
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of weakly compressible turbulence, Sbebslin and Mcatgnary b 
ents of the total pressure (thermal plus magnetic) refax away at the 
Their studies differ from the pres e n t ones in that we are concerne d i 
for which the fluctuations in thermal pr ess ur e and mass density an 
portional. 


In the present calculations, the mass density and pre 

lived or readily discemable trends in spectral oorrefataon, the 
the magnetic field. We can understand this through the fPHy 

(1) In MHD flows, t '» mass density changes in resp o nse ; so total pn 
contributions freon both magnetic and thermal pressures, and 

(2) In compressible flows, the mass density and thermal premier fia 
highly correlated due to die influence of e ntr opy 


The mass density and magnetic press ure did. however, show a net yeiw v- 

eraged) correlation during the evolution of the system at wera g t- Hath anaba Jf = 3-S. 
At such high Mach numbers, compression of the ma gnetic field and the find at 3he *&oafa 
most likely dominates the average correlation (owr the grid) wd roadd he we&cxun so pw 
a net positive corrdaticn, 
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IV- Summary 

The Orszag-Tang vortex system provides a spatially periodic model of a magnetofluid 
undergoing transition to turbulence. The model demonstrates the emergence of small 
scale structure through magnetic reconnection and current sheet formation. Over time, 
dynamic alignment of the momentum and magnetic fields occurs. The Orszag-Tang sys- 
tem also represents an easily reproducible set of initial conditions for comparisons among 
different numerical models of MKD. The present calculations and those described in I 
extend the set of numerical results to compressible MHD flows, including those which are 
moderately supersonic. We have found that the presence of embedded supersonic flows 
and the existence of local force imbalances at moderately high 0 can cause fast shocks and 
broad compressions. The evolution and structure of the “supersonic” Orszag-Tang system, 
therefore, differ significantly from those of its subsonic and incompressible counterparts. 
The same phenomena will occur in any turbulent magnetofluid with embedded supersonic 
regions and local properties similar to the Orszag-Tang system. As a consequence, we 
predict that many structural, global, and spectral properties of such a medium will differ 
substantially from those of subsonic and incompressible turbulent MHD flows. Here we 
summarize those properties. 

In the Orszag-Tang vortex system, the presence of moderately supersonic regions 
{.M s: 2) can result in the emergence of shocks in the region outside i the magnetic X- 
point. Accompanying the shocks are prominent current sheets, which dissipate magnetic 
energy. Such current sheets were not present in our subsonic and incompressible calcula- 
tions. In addition, vortex structures reside in regions of changing shock curvature. The 
shocks occur when the forces due to thermal and magnetic pressure do not balance the 
distribution of momentum in the initial conditions. The excess centrifugal force pulls fluid 
away from the central X-point, rarefying the region and delaying significant reconnection 
there. Since the shock current sheets dissipate magnetic energy at early times, the re- 
connection is weaker as the initial supersonic Mach number increases. The intensity of 
the reconnection increases significantly when the shocks (and the dense fluid which they 
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bound) reach the region of the original X-point. The correlation of the momentum and 
magnetic field begins growing rapidly thereafter. Eventually the growth slows considerably 
so that the correlation reaches a “plateau,” which decreases with increasing Mach number. 
The present model does not indicate with certainty whether reconnection always is en- 
hanced and accelerated by compressive fluid motions and inhibited by rarefying motions. 
However, both effects occur in the present calculations. 

The compressed regions bounded by shocks dominate autocorrelation spectra at early 
times, including that of the solenoidal velocity component. We have not yet explained why 
the latter occurs, other than to infer that the solenoidal and nonsolenoidal components are 
closely coupled and to note that solenoidality of the flow does not imply incompressibility of 
the medium. As opposed to the subsonic and incompressible cases, the normalized spectral 
correlation between the velocity and magnetic fields is near zero over the middle and 
upper portions of the wavenumber domain. The correlation between thermal and magnetic 
pressure is negative over a wide wavenumber range at early times for all calculations with 
M > 1.0. This is consistent with the presence of broad, long-lived regions of compression. 
At later times, these regions are disrupted and the spectral anticorrelation of thermal and 
magnetic pressures disappears. 

To complete our investigation of the compressible Orszag-Tang model, we axe currently 
investigating the effect of varying 0 on the above conclusions. The present and previous 
simulations set 0 > 10/3, so that the magnetic field has not played as dominant a role as 
it might in some space or laboratory plasmas. Our next step in the study of compressible 
MHD turbulence will be to investigate the decay of magnetofluids from random initial 
conditions with and without embedded supersonic flows. 
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APPENDIX A 

Pressure Imbalance in the Supersonic Orszag-Tang Vortex System 


The present model of supersonic regions within the Orszag-Tang vortex system uses 
a fluctuating pressure pi = 6p° which is consistent with Eqs.(6) for an incompressible 
magnetofluid of initial mass density p° = 1. The corresponding velocity and magnetic 
fields, denoted by Vi and Bj, respectively, then satisfy the elliptic equation: 

V . (v, ■ VvO = i [ - i V 2 (p, + |B, I 2 ) + V • (B, • VB, )] . (A.1) 

Equation (A.l) represents the force balance between the motions of the magnetofluid and 
the confining thermal and magnetic pressures. The flat background pressure, p°, and the 
characteristic Mach number, M, must be consistent with p° and |v°| = yjvf), through 
the equation 


7 P° 

For the subsonic cases with p° = 1, Eqs.(A.l), (6), and (9b) represent a consistent speci- 
fication of the initial conditions. As stated in Section II.B, the constraint on the average 
pressure p° limits its minimum value to max{— pi(z, y)} + e, where e > 0, so that the local 
Mach number will be defined everywhere. For 7 = 5/3, Fig. 1 of Dahlburg and Picone 
[1] shows that p° is > 2.6. Equation (9b) then limits the characteristic Mach number 
to M < 0.679. Even when p° is sufficiently close to 2.6 for supersonic regions to exist 
within the initial flow field, we do not observe shocks, and the quantitative behavior of 
such systems is quite close to that of the M = 0.6 case presented in I. 

Equation (9b) shows that one could increase the characteristic Mach number, M, by 
increasing the amplitude of the initial velocity field, increasing the initial mass density 
or lowering p°. If the initial velocity in Eq.(6) were increased by a constant factor, we 
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would also increase the magnetic field amplitude by the same factor in order to maintain 
the initial normalized cross helicity in Eq.(14) at 0.5. Given this constraint, consider each 
alternative separately: 

(1) Scale the velocity field up by a constant factor. Set B' = aBi and v' = avi with 
a > 1. Then by Eq.(A.l), />' = 1, and p' = a 2 pi => p 0 ' > 2.6a 2 . By Eq.(9b), we have 
M < 0.679. 

(2) Scale the mass density up by a constant factor, so that p' = a with a > 1. Then by 
Eq.(A.l), a solution is v' = Vi Jy/a, B' = Bj, and p' = p\ => p 0> > 2.6. By Eq.(9b), we 
have M < 0.679. Note that this solution also fails to maintain the initial normalized cross 
helicity in Eq.(14) at 0.5. 

(3) Scale the initial pressure down by a constant factor, so that p' = ap\ and p 0 ' > 2.6a, 
where a < 1. Equation (A.l) gives the solution v' = y/avi, B' = i/aBj , and p' = 1. 
Again Eq.(9b) shows that M < 0.679. 

This shows that one cannot increase M to supersonic values while maintaining the 
above constraints on v and B and satisfying Eq.(A.l). We can achieve supersonic values 
of M only by violating the balance of forces and pressures implied by Eq.(A.l). Hence, 
our approach is to increase p° above 2.78 and, at the same time, to allow 6p°, v, and B 
to satisfy Eq.(A.l) with p° set equal to one. Since the fluctuating pressure will not balance 
the other local forces, finite magnetoacoustic waves or shock waves will emerge in the flow. 
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Fig. 1 — Contours of constant local Mach number Jl at time t = 0. (a) Characteristic Mach number M = 0.679 with 
balanced forces throughout grid. Minimum value of the initial pressure is 1 x 10' 4 . Minimum and maximum values of 
Jl are 0.0 and 2. 185. (b) Characteristic Mach number M = 1.0 with unbalanced forces in the grid. Minimum and 
maximum values of Jl are 0.0 and 1 .620. 
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Pig 5 — Plots of the generalized stream function for M = 1.0 and S = 50 giving dtreeboo of vdocy ((» »d (O, «d 
magnetic field ((b) and (d)) at respective limes t = 3.0 and 7.0. 

























Fig. 9 — Evolution of global averages for M = 1.0 and 
S = 50 (A) S =100 (O), and S = 200 (□). (a) 
Compressible correlation coefficient a c , (b) Kinetic 
energy, (c) Magnetic energy, (d) Kinetic enstrophy, 
{a 2 }. (e) Magnetic enstrophy, (j 2 ). 









Fig. 10 — (Continued) Autocorrelation spectra (log of magnitude squared) vs. log of squared wavenumber for (a) Mass 
density, t - 1.0, (b) Mass density, t = 2.0, (c) Mass density, t = 4.0, (d) Velocity, t - 2.0, (e) Nonsolenoidal component 
of the velocity, t — 2.0, (0 Solenoidal component of the velocity, t = 2.0, (g) magnetic field, t = 2.0. 





